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ABSTRACT 

A methodology is presented for using the Volterra- 
Wiener theory of nonlinear systems in aeroservoelastic 
(ASE) analyses and design. The theory is applied to the 
development of nonlinear aerodynamic response models 
that can be defined in state-space form and are, therefore, 
appropriate for use in modem control theory. The theory 
relies on the identification of nonlinear kernels that can be 
used to predict the response of a nonlinear system due to 
an arbitrary input. A numerical kernel identification 
technique, based on unit impulse responses, is presented 
and applied to a simple bilinear, single-input-single- 
output (SISO) system. The linear kernel (unit impulse 
response) and the nonlinear second-order kernel of the 
system are numerically-identified and compared with the 
exact, analytically-defined linear and second-order kernels. 
This kernel identification technique is then applied to the 
CAP-TSD (Computational Acroclasticity Program- 
Transonic Small Disturbance) code for identification of 
the linear and second-order kernels of a NACA64A010 
rectangular wing undergoing pitch at M=0.5, M=0.85 
(transonic), and M=0.93 (transonic). Results presented 
demonstrate the feasibility of this approach for use with 
nonlinear, unsteady aerodynamic responses. 

INTRODUCTION 

The subject of nonlinear unsteady aerodynamics is one 
of great interest in the aerospace community. The interest 
is due to the fact that nonlinear unsteady aerodynamic 
behavior can have a significant effect on the performance 
and stability of a flight vehicle. It is therefore very 
important to be able to predict and understand nonlinear 
unsteady aerodynamic behavior. 

Today's most powerful and sophisticated tools for 
predicting nonlinear unsteady aerodynamics are being 
developed in the field of computational fluid dynamics 
(CFD)l. The nature and detail of the nonlinear fluid flow 
that is predicted by a particular flow solver depends on the 
governing equations that are discretized in the solver. The 
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order of the governing flow equations can vary from the 
transonic small disturbance (TSD) level to the full Navier- 
Stokes equations. As CFD methods improve our ability 
to predict nonlinear unsteady flows, it is a natural and 
important step to investigate methods for controlling 
these flows in order to improve the performance and/or 
stability of a flight vehicle at flight conditions where 
nonlinear unsteady aerodynamic effects are significant 
Modem aeroservoelastic (ASE) analysis tools, such as 
ISAC^ (Interaction of Structures, Aerodynamics, and 
Controls) and ADAM*^ (Analog and Digital 
Aeroservoelastic Method), are routinely used for predicting 
the interaction between the structural system, the 
aerodynamic system, and the control system of a flexible 
aircraft so that control laws that account for and take 
advantage of this flexibility can be designed. The goal of 
the control law design may be for flutter suppression 
(stability)^ and/or for load alleviation (performance), but, 
in either case, the control system design has been limited 
to linear aerodynamic responses. This limitation inhibits 
the design and analysis of control systems that can 
account for flow nonlinearities such as shocks, boundary 
layer effects, and separated flows. Although nonlinear 
aerodynamics are eventually incorporated into the control 
system design via wind-tunnel studies or semi-empirical 
simulations, there is a very real need for modeling 
nonlinear aerodynamic behavior, such as that predicted by 
CFD codes, early in the design phase for use in ASE 
analysis methods. Although some work has been done in 
directly incorporating simple control laws into CFD 
codes^*\ these approaches do not generate a general 
model of the nonlinear aerodynamic response. Instead, the 
control law gains have to be varied in a trial-and-error 
manner as flight conditions are varied. Since linear 
aerodynamic response is modeled as a linear system using 
rational function approximations in modern ASE codes, 
the purpose of this research is to investigate the 
feasibility of modeling nonlinear aerodynamic response as 
a nonlinear system, in particular, as a Volterra-Wiener 
nonlinear system. 

The Volterra theory was developed by Volterra in 
1930^. The theory is based on functionals, or functions 
of other functions, and subsequently became a 
generalization of the linear convolution integral approach 
that is applied to linear time-invariant (LTI) systems. 
The theory formulates the response of a nonlinear, time- 
invariant system as an infinite sum of multidimensional 
convolution integrals of increasing order, with the first 
term in the scries being the standard linear convolution 
integral. Each multidimensional convolution integral in 



the series has an associated kernel. The first-order kernel 
is simply the linear unit impulse response of the system 
and the higher-order kernels are measures of nonlinearity 
of the system response. This infinite sum of 
multidimensional convolution integrals is known as the 
Voltenra Series and it is well defined in both the time and 
frequency domains. 

The Volterra theory has been applied primarily to 
nonlinear electrical and electronic systems. Wiener^ 
contributed significantly to the Volterra theory and, as a 
result, the theory is currently referred to as the Volterra- 
Wiener theory of nonlinear systems. References 9 and 10 
developed a kernel identification technique based on auto- 
and cross-correlation functions that can be applied to 
nonlinear, time-varying systems. The textbooks by 
Rugh* 1 and Schetzen^ and the work by Boyd, Chua, and 
Desoer^ and several others 14-18 ^ excellent, detailed 
descriptions of the Volterra-Wiener theory and are highly 
recommended to the interested reader. 

Application of nonlinear system theories, including 
the Volterra-Wiener theory, to the problem of modeling 
nonlinear unsteady aerodynamic responses has not been 
extensive. Ueda and Dowell’s 1 ^ application of the 
concept of describing functions to unsteady transonic 
aerodynamic responses is one approach. The work by 
Tobak and Pearson^® involved the application of 
Volterra’s concept of functionals to indicial aerodynamic 
responses for the analytical derivation and experimental 
determination of nonlinear stability derivatives. The work 
by Jenkins^l is also an investigation into the 
determination of nonlinear aerodynamic indicial responses 
and nonlinear stability derivatives. Stalford, Bauman, 
Garrett, and Herdman22 successfully developed Volterra 
models for simulating the nonlinear behavior of a 
simplified nonlinear stall/post-stall aircraft model and the 
behavior of a simplified model of wing rock. 

In Ref. 22, the nonlinear aerodynamic response is 
analytically defined a priori so that derivation of the 
Volterra kernels is a straightforward procedure. In general, 
the nonlinear response of a given configuration at a given 
flight condition will not be known. The output from a 
CFD code provides information regarding the nonlinear 
aerodynamic response of the configuration at a given 
flight condition to a given input, but a limited amount of 
information can be inferred regarding the nonlinear 
aerodynamic response of the configuration to an arbitrary 
input. Prediction of the nonlinear aerodynamic response 
of a configuration to an arbitrary input requires 
identification of the nonlinear kernels of the Volterra 
Scries for the particular configuration. 

The problem of Volterra kernel identification has been 
addressed by Rugh 11 , Clancy and Rugh 2 3 for discrete 
systems, Schclzen^ 4 , and more recently Boyd, Tang, and 
Chua^S. There are several ways of identifying the 
Volterra kernels, both in the time and frequency domains. 
The methods can be applied to continuous or discrete 
systems^, such as CFD models. Recently, Tromp and 
Jenkins^ applied a Laplace domain scheme and 
aerodynamic indicial responses for the identification of the 
Volterra kernels of a two-dimensional airfoil undergoing 


pitching motions using a Navier-Stokes flow solver. The 
first order, or linear, kernel was identified. The second- 
order kernel was identified for a sample problem and the 
method of Boyd, Tang, and Chua^S W as suggested for 
identification of the second-order nonlinear kernel of the 
airfoil response. 

An important characteristic of the Volterra-Wiener 
theory of nonlinear systems is that a bilinear state-space 
system can be realized once the nonlinear kernels of 
interest have been identified27»28 This bilinear state- 
space system can be used as a nonlinear aerodynamic 
model for aeroservoelastic analysis and design. Standard 
control theory techniques or the theory of bilinear optimal 
control can then be used for designing control systems 
that account for nonlinearities in the aerodynamic 
response. 

The objective of the present research is to investigate 
the application of a time-domain identification technique 
to the three-dimensional CAP-TSD (Computational 
Aeroelasticity Program-Transonic Small Disturbance) 
code29 for identification of the nonlinear, second-order 
kernel of a NACA64A010 rectangular wing undergoing 
pitch at transonic Mach numbers. The paper begins with 
a brief description of the Volterra-Wiener theory of 
nonlinear systems followed by the description of a kernel 
identification technique based on unit impulse responses. 
The kernel identification technique is applied to a simple 
bilinear state-space system to provide insight into the 
application of the technique and the nature of a nonlinear 
kernel. The CAP-TSD code and the computational model 
of the rectangular wing are then described and, finally, 
some results for the wing are presented and discussed. 

VOLTERRA-WIENER THEORY 

The basic premise of the Volterra-Wiener theory of 
nonlinear systems is that any nonlinear system can be 
modeled as an infinite sum of multidimensional 
convolution integrals of increasing order. This infinite 
sum is known as the Volterra Series and it has the form 

t 

y(l) = J h x (t- x) u(t) dr + 

0 
t t 

jj h^t - x,,t - T 2 ) u(x,) U(x 2 ) dx, dx 2 + ... 

00 

L t 

+ j-J hm(t - - In) u(T,)...u(T n ) dt,...dX n + ... (1) 

0 0 

where y(t) is the response of the nonlinear system to u(t), 
an arbitrary input; hi is the first-order kernel or the linear 
unit impulse response; h2 s is the second-order kernel, and 
h ns is the nth order kernel. It is assumed that; 

1). the kernels, input function, and subsequendy, the 
output function are real-valued functions defined for 
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tj e (-oo,+oo), i=l,...,n,... 

2) . the system is causal so thal h ns (X] ,...,x n ) = 0 if any 
Xj<0 

3) . the system is time invariant 

Inspection of Equation (1) reveals some very 
interesting and characteristic features of the Volterra 
series. If the kernels of order two and above arc zero, then 
the response of the system is linear and is completely 
described by the unit impulse response hi(t), and the first 
order convolution integral. The assumption underlying 
the first order, or linear, convolution integral is that the 
response of the system at a given time, l, is independent 
of the response of the system at a previous lime. This is 
why convolution of a single unit impulse response with 
an arbitrary input is valid for predicting the response of a 
linear system. The higher order kernels h nSj are the 
responses of the nonlinear system to multiple unit 
impulses, with the number of impulses applied to the 
system equal to the order of the kernel of interest : c.g., 
h2s is the response of the nonlinear system to two unit 

impulses applied at two varying points in time, x\ and 
X2. The mathematical definition follows directly for the 
nth order kernels although visualization of these functions 
can become difficult for orders greater than three. The 
nonlinear kernels are measures of the relative influence of 
a previous input on the current response, which is a 
measure of nonlinearity. This temporal measure of 
nonlinearity is referred to as memory and, as a result, 
Volterra systems are sometimes referred to as nonlinear 
sytems with memory. 

The 's' in the kernel names stands for 'symmetric' 
since h2 s (xi»X2) = h 2s(X2Ji). Although, depending on 
the domain of integration that is chosen, the kernels can 
be defined in 'triangular' or 'regular' form, any kernel can 
be symmetrized without affecting the input/output 
relation. This is done by realizing that 

hsym(xi>-»»x n )=(l/n!) £ h(Xn(i Xn( n )) (2) 

where the indicated summation is over all n! permutations 
of the integers 1 through n. For the present study, only 
symmetric kernels will be investigated since these are 
mathematically easier to interpret and intuitively easier to 
visualize. For the interested reader, details regarding this 
issue can be found in Refs. 1 1 and 12. 

One approach for obtaining Volterra series 
representations of physical systems is to assume that the 
system is a 'weakly' nonlinear system. A system that is 
weakly nonlinear is a system that is well defined by the 
first few kernels of the Volterra scries so that the 
magnitudes of the kernels greater than second or third 
order fall off rapidly and are negligible. Boyd, Tang, and 
Chua^5 mention some physical systems that arc 
accurately modeled as weakly nonlinear systems including 
electromechanical and electroacoustic transducers and some 
biological systems. In this initial study, it is assumed 


that the nonlinear aerodynamic system that is synthesized 
from the transonic small-disturbance (TSD) potential 
equation is a weakly nonlinear, second-order system. 
Results are therefore, limited to the identification of the 
second-order kernel, or h2 s . 

It should also be noted that the kernels, linear and 
nonlinear, are input dependent. For example, for a linear 
system, if the response of the system to an arbitrary input 
is desired, the unit impulse response of the system due to 
that particular type of input must first be defined. For a 
single-input-single-output (SISO) system, there is only 
one unit impulse response. For a multiple-input- 
multiple-output (MIMO) system, there are nxm unit 
impulse responses where n is the number of inputs and m 
is the number of outputs. These unit impulse responses 
arc then combined to form the unit impulse matrix. 

Kernel I dentification 

The advantage of the Volterra series approach for 
modeling nonlinear systems is that once the kernels are 
identified, the response of the nonlinear system to an 
arbitrary input can be predicted. The problem of kernel 
identification, therefore, is central to the successful 
generation of an accurate Volterra series representation of 
a nonlinear system. The most obvious approach for 
identifying the kernels is to derive analytical expressions 
for the kernels from the governing nonlinear equations of 
the system of interest^* 31 . Although this approach is 
applicable to any set of nonlinear equations, including the 
nonlinear fluid flow equations such as TSD, Euler, and 
Navier-Stokes equations, it would require some additional 
coding in order to numerically compute the kernels. 
Instead, a kernel identification technique is desired that 
uses the output of a CFD code directly for quick and 
efficient kernel identification. 

Boyd, Tang, and Chua^ describe a frequency-domain 
technique that was successfully applied to the 
experimental identification of the second-order kernel of a 
nonlinear electro-acoustic transducer (speaker) system. 
This and other frequency-domain techniques are available, 
but it is preferable to use a time-domain kernel 
identification technique since unsteady, nonlinear CFD 
analyses are generally performed in the time domain. The 
time-domain method investigated in this study is the 

method of unit impulse responses**’ 12,24 Although 

unit impulse responses are defined for continuous 
systems, Clancy and Rugh^ 3 have shown that an 
equivalent technique using the unit pulse response can be 
used for the kernel identification of discrete nonlinear 
systems. 

In what follows, the kernel identification technique 
using unit impulse responses is derived. The technique is 
then applied to a simple problem in order to illustrate the 
discrete application of the technique and the nature of the 
second -order kernel that is identified. 

A weakly nonlinear, second-order system is described 
by 
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( 5 ) 


y(t) = J h,(t - x) u(t) dx + 

0 

t t 

J J 1^(1 - x, ,t - x 2 ) u(Tj) u(x 2 ) dx, dx 2 (3) 

00 

Inputs consisting of single and double impulse functions 
can be defined as 

uo(0 =8o(t) 

X- 

Ui(t) = 6o(t) + 5o(t - Ti) 

where Ti is a distinct positive number. The responses of 
the system to these two inputs are 

yo(0 = MO + h 2 S (t ,0 

yi(0 = MO + M l - T i) + MCi.O + 

2h 2 s(t.t-Ti) + h 2 s(t-Ti.l-Ti) 

The 2h2s(t. I - Ti) term is a result of the symmetry of the 
kernel since h 2 s (t, t - Tj) = h 2 s (t - Tj ,t). Then 

yi(0 = yo(0 + hi (t - Ti) + 2h2s(l, t - Tj) + 
h2s(t - Tj, t - Ti) 

and noticing that 

yo(t-T 1 ) = h,(t-T 1 ) + h 2 s(l-T l .l-Tj) 

results in 

yi(0 = yo(0 + yo(t - Ti) + 2h 2s (t, t - to 

Solving for the second-order kernel 

h 2s (t,t-T0 = (l/2)(yi(0- yo(0 - yo(l - TO) (4) 

which is the value of the second-order kernel for any value 
of Tj. 

The procedure for computing h 2s is as follows. First, 
yo(t), which is the response of the system to a unit 
impulse response applied at time t, is generated. Then, 
since the system is time invariant, yo is shifted in time to 
a new time t - Tj, which becomes yo(t - TO- The 
response of the system to two unit impulses, one at time 
t and one at time t - Tj, or yi(t), is generated and finally, 
all three responses are substituted into equation (4). The 
second-order kernel, lt2 s is a two-dimensional function of 
time. That is, it is a function of time l and a function of 
time lag Ti so that for every value of T] that is used, a 
new function of lime t is defined. These sets of functions 
of time are referred to as "terms" of the kernel. The first 
term of h2s is defined when Tj = 0, or when both unit 
impulse inputs are applied at the same point in time. 
When Ti=0, equation (4) reduces to 


h 2s (t. 0 = (l/2)(yi(0 - yo(0 - yo(0) 

= 0/2)yi(t) - yo(t) 

The second term of the kernel depends on the next value 
of Ti selected. The number of T i 's, or the number of 
terms, needed to accurately define a second-order kernel 
depends on the nonlinear system under investigation. 

In addition, the linear portion of the nonlinear 
response can be identified when Ti = 0. It is important to 
realize that the linear portion of the nonlinear response is 
not, in general, equivalent to the purely linear response. 
For example, for an aerodynamic system, the linear 
response computed using the linear equations (fiat plate) 
is not identical to the linear portion of the response 
computed using the nonlinear equations (thickness). 

The linear portion of the nonlinear response is defined 
as follows. The response of the system represented by 
equation (3) to 2uo(t) is 

y 2 (t) = 2hi(t) + 4h2 s (t,t) 

Then, solving simultaneously with yo(t) results in 

hi(t)=2y 0 (t) - (l/2)y 2 (t) (6) 

which is the unit impulse response of the linear portion 
of the nonlinear response. 

The equations derived above for h2 s are clear and 
simple measures of nonlinearity. In these equations, 
nonlinearity is measured as a deviation from linear 
behavior. This is evident in that h 2 s is identically zero 
for a linear system due to the principle of linear 
superposition. Therefore, an additional benefit of the 
second-order kernel is that it can be used to measure the 
true linearity of a system that is classified as being linear 
or for establishing boundaries beyond which the 
assumptions of linearity begin to fail. Definitions of 
higher order kernels can be derived in the same way as for 
h2s by applying the appropriate number of unit impulses 
to the system. 

Once h2s is identified, the nonlinear response of the 
weakly nonlinear, second-order system to an arbitrary 
input can be determined. A further advantage of the 
Volterra theory of nonlinear systems is that a bilinear 
state-space system can be realized once the kernels are 
identified, generating a nonlinear aerodynamic state-space 
system that is amenable for use with modem control 
theory. The relationship between the Volterra kernels and 
the bilinear state-space formulation is as follows. 

It is well known that, for a linear system described by 
the following state-space representation 

x = Ax + Bu 

y = Cx (7) 

(where the D matrix, or feedthrough term, has been set to 
zero) the system's unit impulse response is given by 
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h(t) = C (exp(At)) B 


( 8 ) 


If the unit impulse response of a linear system is known, 
then a state-space realization can be obtained using 
standard linear system realization techniques. In a natural 
extension to this concept, a Voltcrra system can be 
modeled by a bilinear state-space system as follows 

x = Ax + Nxu + Bu 

y = Cx (9) 

where the matrix N is a measure of the nonlinearity of the 
system. The system is linear when N=0, as it is reduced 
to equation (7). The relationship between the second-order 
Volterra kernel and the bilinear stale-space formulation is 

h2reg(tl,t2) = C (exp(Ati)) N (cxp(Al2)) B (10) 

where the subscript 'reg' stands for the 'regular' definition 
of the kernel. Equation (2) can then be used to compute 
the symmetric kernel h2 s from h2rcg- 

EXAMPLE PROBLEM 

Assume a system described by the following 
differential equation 

y(0 + a,y(t) + a^t) = b Q u(t) + n 0 y(t)u(l) (11) 

In bilinear state-space form, equation (1 1) can be rewritten 
as equation (9) where 



" 0 1 


0 

0 


" 0 " 

A = 


. N = 



. B = 



1 

V 

cp° 
1 


n o 

0 


_ V 


and C = [ 1 0] 

Values for the A, N, and B matrices were arbitrarily 
chosen as ao~2., aj=3., no =-6., and bo=l. The linear 
response of this system is obtained by setting no = 0. An 
arbitrary input, shown in the inset of figure 1, was defined 
as 

u #rb (t) = 0.0 , t < . 1 and t > 1 . 1 
= t-.l , .1 <t< 1.1 

and applied to discretized versions of the linear (N=0) and 
nonlinear (bilinear) systems in order to examine the 
differences in the responses due to this arbitrary input 
using a time step of 0.01. The responses can be seen in 
figure 1 where it is clear that the effect of the added 
nonlinearity is to reduce the level of the linear response. 

Application of the kernel identification technique to 
continuous systems requires the use of unit impulse 
inputs. For discrete systems, the equivalent of a unit 
impulse input is a unit pulse input defined as 


Up(t) = 1 . , t = tO 
= 0. , t * to 

where tO can be any time step since the system is time 
invariant. Also, the unit impulse response of a 
continuous system is approximately equal to the unit 
pulse response of the discretized version of the same 
system divided by the time step used in the discretization. 
The discrete responses are therefore divided by the time 
step so that comparisons with the continuous, or 
analytical, responses can be made. 

The unit pulse input, Up(t), was applied to the 
discretized linear system and the resultant unit pulse 
response was compared to the analytical (or continuous) 
unit impulse response of the system from equation (8). 
This comparison is presented in figure 2, showing 
excellent agreement. The slight difference between the 
two responses is a result of the time step (DT=0.01) that 
was used. A smaller time step improves the accuracy of 
the unit pulse response. 

The unit pulse input was then modified to include two 
unit pulses for identification of the second-order kernel so 
that 

up(t)= 1. , t = tO and t = tl 
= 0. , t * tO and t * tl 

where the value of tl is a time step such that tl-tO=T|. 
The tO was held fixed while the tl was varied in order to 
vary Ti- The value of Ti was varied in increments of ten 
time steps, or 0.1 seconds. The first term of the second- 
order kernel corresponds to Ti=0.0; the second term of the 
second-order kernel corresponds to Tj=10 time steps ; the 
third term corresponds to Ti=20 time steps and so on. A 
total of twenty terms were generated and the resulting 
second-order kernel is shown in figure 3 as a function of 
time t and time lag Ti. In figure 3, TB is the second 
term of the kernel (T]=10 time steps), TC is the third 
term, TD is the fourth term, and the first term, TA, is not 
visible because it was zero. As can be seen, the second- 
order kernel is a two-dimensional function that varies 
smoothly with time t and time lag Tj. It can also be 
seen that the second-order nonlinear response of the 
system at each time lag T] exhibits a rapid initial growth, 
reaches a maximum response and then begins to dissipate. 
The fact that the second-order kernel is negative is 
consistent with the result in figure 1 where it was shown 
that the added nonlinear terms reduced the response of the 
system from the purely linear response. Inspection of a 
nonlinear kernel, therefore, can provide a significant 
amount of information regarding the behavior of a 
nonlinear system in terms of amplitude and time lag 
variations. This is important information for the design 
of effective control systems. 

The exact, analytical second-order kernel of the 
bilinear system was computed using equation (10). The 
analytical second-order kernel is not a symmetric kernel 
and so it must be symmetrized before it can be compared 
with the numerically-identified symmetric second-order 
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kernel. Symmetrization of ihc analytical kernel is 
performed by using equation (2). The numerically- 
identified, or discrete, kernel is divided by the square of the 
lime step so that it can be compared with the analytical, 
or continuous, kernel. The symmetrized analytical 
second-order kernel and the symmetric, numerically- 
identified second-order kernel are shown in figure 4 for 
two time lags, TD (30 time steps) and TG (60 time 
steps). The comparison is excellent with slight 
differences occuring around the regions of maximum 
response. Again, improved accuracy can be achieved by 
using a smaller time step in the numerical identification 
technique. The important point to be made here, 
however, is that the use of discrete, unit pulses can be 
used to identify the discrete, second-order Vollcrra kernel 
of a discrete nonlinear system. 

COMPUTATIONAL PROCEDURES 

The CAP-TSD program is a finite-difference program 
which solves the general-frequency modified transonic 
small-disturbance (TSD) potential equation. The TSD 
potential equation is defined by 

m 2 (<t>, + 2<t>, ) t = [(1 -m! )<K + F<t>, + C<t> y 2 | x + 

(<|> y + HtMy ) y + (<|> Z \ (12) 

where Moo is the freestream Mach number, <)> is the 
disturbance velocity potential, and the subscripts of <|) 
represent partial derivatives. 

Several choices are available for the coefficients F, G, 
and H depending upon the assumptions used in deriving 
the TSD equation. For transonic applications, the 
coefficients are herein defined as 

F =-j(Y+1)M» 

G = j(Y-3)m! 

H=-(y-1)m2 (13) 

where y is the ratio of specific heats. The linear potential 
equation is recovered by simply setting F, G, and H equal 
to zero. 

Equation (12) is solved within CAP-TSD by a lime- 
accurate approximate factorization (AF) algorithm 
developed by Batina^, in Refs. 32 and 33, the AF 
algorithm was shown to be efficient for application to 
steady or unsteady transonic flow problems. Several 
algorithm modifications have been made which improve 
the stability of the AF algorithm and the accuracy of the 
results^ 4 . One of these improvements is the option to 
include vorticity and entropy corrections^ for improved 
shock modeling. This option, however, was not used for 
the analyses presented in this paper. 


The CAP-TSD program can treat configurations with 
combinations of lifting surfaces and bodies including 
canard, wing, tail, control surfaces, tip launchers, pylons, 
fuselage, stores, and nacelles. The code was recently 
applied to the Active Flexible Wing (AFW) wind-tunnel 
model, which included modeling of the fuselage and tip 
stores, for prediction of the model's transonic aeroelastic 

behavior^. 

The code has an exponential pulse capability that can 
be used for generating unsteady aerodynamic pitch, 
plunge, and modal responses. The pulse and pulse rate are 
defined as 

p(t) = 8 0 exp(-w (t - tc) 2 ) (14) 

p(t) = -2w (t - tc) p(t) (15) 

where t and tc are in terms of nondimensional time steps. 
For pitching motions, the angle-of-attack input function 
is defined using equation (14) and the rate of change of 
angle of attack, or angle-of-attack rate is defined using 
equation (15). These functions of time then become part 
of the downwash equation which, for simple pitching 
motions is 

f'tx.t)* = -^2- - a (l) - a (t) (x - x piich )/tL ( 1 6) 


where the plus and minus signs indicate upper and lower 
surfaces of the airfoil. The first term in equation (16) is 
the airfoil geometry slopes, followed by the angle of 
attack, and by the angle-of-attack rate multiplied by the 
pitch distance where Xpi tc h * s ^e pitch ax * s - Th e 
downwash provides the boundary condition defined at the 
z=0 plane required to complete the solution of the TSD 
equations. For the linear aerodynamic solution, or flat 
plate solution, the airfoil geometry slopes are zero and the 
downwash becomes a function of angle of attack and 
angle-of-attack rate only. 

The computational grid of the NACA64A010 
rectangular wing is dimensioned 137 by 40 by 84 grid 
points in the x-, y-, and z-directions respectively. The 
wing has an aspect ratio of three but the computational 
domain covers only one semi-span due to flow symmetry. 

RESULTS AND DISCUSSION 

Before the kernel identification technique was applied 
to an aerodynamic system, the appropriate excitation 
input had to be defined. It was obvious that the excitation 
input had to be a perturbation of the downwash, equation 
(16). The excitation input could, therefore, be composed 
of one or a combination of the parameters that define the 
downwash. In addition, the excitation input had to be of 
an "impulsive" nature so that the theoretical assumptions 
used in the derivation of the kernel identification technique 
were maintained. 

The requirement that the excitation input be 
"impulsive" disallowed the use of the exponential pulse 
capability defined in the CAP-TSD code, eqs. (14) and 
(15). The aerodynamic response induced by an 
exponential pulse input is a response of the aerodynamic 
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system to a smoothly-varying function of angle of attack 
and a smoothly-varying function of angle-of-attack rate. 
This response cannot be referred to as the unit pulse 
response of the system and therefore cannot be used in the 
linear convolution integral to predict the linear response 
of the system or in the identification of the nonlinear 
kernels. The correct excitation input should, therefore, be 
of a unit magnitude and should be applied at only one 
time step, as was presented for the example problem 
earlier. 

Based on this reasoning, the unit pulse inputs 
available are then 

a (t) = 0.01745 rads (1 deg) , t = tO 
= 0.0 , t * tO 

with fit (t) = 0.0 everywhere 
or, 

Ct (t) = 1 rad/sec , t = tO 
= 0.0 , t * tO 

with a (t) = 0.0 everywhere 

or a combination of both angle of attack and anglc-of- 
attack rate inputs. Also, due to the impulsive nature of 
the unit pulse inputs, a very small time step of 
DT=0.0001 had to be used in order to obtain smooth 
aerodynamic responses. All dynamic responses, linear and 
nonlinear, were 500 time steps in length. It should be 
noted that this choice of time step and number of time 
steps results in only .05 chords of travel. This very small 
time sample is dominated by the high frequency content 
of the response so that the results that follow are limited 
to high frequency responses. 

The application of the unit angle-of-attack pulse input 
resulted in very small aerodynamic responses for both the 
linear (flat plate) and the nonlinear (thickness) 
aerodynamic solutions using CAP-TSD. An attempt to 
identify the nonlinear kernel using these responses 
resulted in numerical noise and was not possible. The 
second type of input, the unit angle-of-attack rate pulse 
input, provided sufficient excitation of the nonlinear 
equations so that a nonlinear kernel could be identified, as 
will be seen. The combined unit angle-of-attack and unit 
angle-of-attack rale inputs resulted in responses that were 
only marginally different from the unit angle-of-attack rate 
responses. As a result of these preliminary 
investigations, all subsequent analyses arc based on unit 
angle-of-attack rate pulse inputs only. Also, only lift 
coefficient responses to pitching motions about the 
quarter-chord location are investigated in this study, 
although the techniques presented arc applicable to any 
force coefficient response including generalized 
aerodynamic forces. 

For the results that follow, PI is the unit pulse 
response at time tO; P2 is the unit pulse response at time 
t=t0+Ti; P12 is the response due to a unit pulse at time 


t=t0 and another unit pulse at time t=t0+Ti; and Pll is 
the response due to two unit pulses at time t=t0. The 
unit angle-of-attack rate inputs were applied at the 60th, 
1 10th, and 160th time steps which translates to time lags 
equal to 0.0, 50.0, and 100.0 time steps. The choice of 
these time steps was arbitrary and the choice of the 60th 
time step as the first response, or as the PI response, was 
done in order to avoid any numerical transients that might 
occur when the nonlinear aerodynamic analyses are 
initiated from steady-state solutions. 

Linear Kernel Identification 

The PI linear (fiat plate) lift-coefficient response 
shown in figure 5 has the characteristics that are typical of 
a unit pulse response. That is, the initial part of the 
response is impulsive and the latter part of the response is 
a damped transient. It should be restated that the time 
step at which the input is applied is not important since 
the linear system is time invariant 

The linear convolution integral, 

'i 

y(0 = J h(l - x) u(x) dx 

0 

was used to verify that PI (figure 5) was indeed a unit 
pulse response. This was done by first generating a linear 
(flat plate) lift-coefficient response to an arbitrary input 
using the CAP-TSD exponential pulse capability, eqs. 
(14) and (15) at M=0.85. The pulse was defined with a 
w=90,000.,8o = 0.009 rads (.5 deg), and centered at 
tc=250 lime steps. The response to this exponential 
pulse, which will be refered to as the linear arbitrary 
response, is shown in figure 6 along with the angle of 
attack (inset) and the corresponding angle-of-attack rate 
functions generated by eqs. (14) and (15). It should be 
noted that this response is of a relatively high-frequency 
which, as can be seen in figure 6, is influenced primarily 
by the angle-of-attack rate input. 

The linear convolution integral was then evaluated 
using the angle-of-attack rate function, presented in figure 
6, as the input, u, and the PI response shown in figure 5 
as the unit impulse response, h. A comparison of the 
linear arbitrary response shown in figure 6 and the 
response computed from the linear convolution of PI and 
the angle-of-attack rate function is presented in figure 7. 
The excellent agreement between the two responses 
verifies the use of the unit angle-of-attack rate pulse input 
for generating unit pulse responses that can be used for 
predicting high-frequency arbitrary responses. However, 
accurate prediction of low-frequency responses, where the 
angle of attack input becomes significant, needs further 
investigation. 

Nonlinear Kernel Id entification 

The nonlinear aerodynamic responses were computed 
about steady-state converged solutions of the 
NACA64A010 rectangular wing at zero degrees angle of 
attack. These steady-state solutions consisted of about 
2500-5000 time steps using a time step DT=0.01. The 
gas medium used in the nonlinear analyses was air which 
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corresponds to a 7=1.4. The assumption of time 
invariance was verified for the nonlinear aerodynamic 
responses. 

At M=0.85, in the steady solution, a strong shock is 
present at the 60% chord location so that a level of 
nonlinear response should be noticed when identification 
of the second-order kernel is performed. The nonlinear 
PI, P2, P12 and Pll responses of the NACA64A010 
wing at M=0.85 are shown in figure 8. The linear 
portion of the nonlinear response, or h] (eq. (6)), is 
computed first as 

hi = 2(P11) - .5(P1) 

The linear unit pulse response from the linear (flat plate) 
solution is ’h 1 and should not be confused with 'hi'. 
Figure 9 is a comparison of these two unit pulse 
responses where the difference is small but noticeable, in 
terms of an associated time lag for the hi . 

The first term of the second-order kernel was computed 
using 

H2S1 = .5(P11-P1-P1) = .5(P1 1) - PI 

which corresponds to a time lag Tl=0.0. The nonlinear 
P2A and PI 2 A responses were computed using a time lag 
Ti=50.0 and the nonlinear P2B and P12B responses were 
computed using a time lag Ti=100.0. The P2A, P12A, 
P2B, and P12B responses were used for computing the 
second and third terms of the second-order kernel, 
respectively, as 

H2S2 = .5(P12A - P2A - PI) 

H2S3 = ,5(P12B -P2B - PI) 

Notice that computation of additional terms of the second- 
order kernel requires only the generation of the appropriate 
PI 2 response since the PI response needs to be computed 
only once and the P2 response is just ihe PI response 
shifted in time. 

The three terms of the M=0.85 second-order kernel, 
H2S1, H2S2, and H2S3, are shown in figure 10. It can 
be seen that the identified second-order kernel, although 
noisy at the smaller magnitudes, docs exhibit a particular 
shape which is not numerical noise. The small values of 
the M=0.85 second-order kernel arc indications that the 
nonlinearities at this condition arc small. This is also 
consistent with the result presented for the example 
problem where the second-order kernel identified for that 
case was on the order of l.c-06. Also, the approach to 
zero values of the second-order kernel as time lag is 
increased, or as more terms are computed, is as expected 
since the second-order kernel is a finite and bounded 
function of time. The M=0.85 second-order kernel 
therefore exhibits behavior that is characteristic of a 
second-order kernel. 

The same three terms of the second-order kernels at 
M=0.5 and M=0.93 were also identified. At M=0.93, the 
shock was located at the trailing edge. Figure 1 1 is a 
comparison of the nonlinear unit pulse responses of the 


NACA64A010 rectangular wing at M=0.5, 0.85, and 
0.93. It is noticed that the amplitude and frequency of 
these unit pulse responses decreases as Mach number is 
increased. This, of course, is configuration and motion 
dependent so that for a plunging motion, for example, 
this same trend may not occur. Comparison of the 
second-order kernels for all three Mach numbers reveals a 
suprising and interesting result. Figure 12 presents the 
three terms of the M=0.5 second-order kernel and figure 13 
presents the three terms of the M=0.93 second-order 
kernel. At first glance, and comparing with the three 
terms of the M=0.85 second-order kernel (figure 10), it 
appears that the magnitude of the second-order kernels 
decreases with increasing Mach number. This is the 
reverse trend that is expected since the second-order kernel 
should increase with increasing levels of nonlinearity, or, 
in this case, with an increase in Mach number. 

This preliminary comparison is not complete. An 
appropriate comparison of the second-order kernels of the 
three Mach numbers should account for the differences in 
the magnitude of the unit pulse responses at each Mach 
number. If, for example, the maximum absolute value of 
each of the three terms of the second-order kernels is 
divided by the maximum value of the unit pulse response 
of each corresponding Mach number, the resulting values 
would be better indicators of the relative magnitude of the 
nonlinear effects. A bar chart of these values, referred to 
as the relative nonlinearity, for the three Mach numbers is 
presented in figure 14 for all three terms of the second- 
order kernels. Figure 14 does indeed show the growth of 
relative nonlinearity as Mach number is increased for all 
three terms of the second-order kernels and, in particular, 
the sudden increase in the first term at M=0.93. 

A qualitative interpretation of the second-order kernels 
identified at M=0.5, M=0.85, and M=0.93 would imply 
that a nonlinear response at these Mach numbers is 
dominated by an in-phase increase in magnitude (the first 
term) and subsequent, but smaller, time-lagged variations 
in the response. This is, in fact, a reasonable 
interpretation as can be seen in figure 15, which is a 
comparison of the linear (flat plate) aerodynamic and 
nonlinear (thickness) aerodynamic responses at M=0.85 
due to the same exponential pulse function described in 
figure 6. The nonlinear response is larger in magnitude at 
the regions of maximum response and exhibits some 
phase difference with respect to the linear response. 

The results presented thus far are encouraging and 
support the feasibility and applicability of the Volterra 
series approach for the modeling of nonlinear unsteady 
aerodynamic responses. Additional work is needed for 
determining the minimum number of terms of the second- 
order kernel that are required for accurate prediction of 
responses to arbitrary inputs, in applying the unit pulse 
response technique to lower frequency responses where the 
effect of angle of attack inputs will be more significant, 
and in evaluating the third order kernels for verification of 
the assumption of a weakly nonlinear system. 

CONCLUSIONS 

The Volterra- Wiener theory of nonlinear systems was 
briefly described and presented as a method for modeling 
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nonlinear aerodynamic responses for use in 
aeroservoelastic (ASE) analysis and design. Successful 
application of the theory is contingent upon the sucessful 
identification of the nonlinear kernels. Although several 
methods for identifying the kernels exist, the method 
investigated in this study was a time-domain technique, 
based on unit impulse responses. The technique was 
described and applied to a simple bilinear (nonlinear) 
system for illustrative purposes and to verify the 
application of the technique to discrete nonlinear systems. 
The first order, or linear, kernel and the second-order 
nonlinear kernel of the simple system were accurately 
identified using the technique. Interpretations of the 
second-order kernel were also presented. 

The kernel identification technique was then applied to 
a NACA64A010 rectangular wing using the CAP-TSD 
code. Application of the kernel identification technique to 
the aerodynamic model began with the identification of 
the first order kernel due to the linear aerodynamic 
response (flat plate solution) of the wing in pitch about 
the quarter-chord location. It was shown that a unit 
impulse response, or, for discrete systems, a unit pulse 
response, can be accurately computed by a unit excitation 
of rate of angle of attack in the downwash function for 
predicting linear high-frequency responses. Additional 
work needs to be performed for prediction of low 
frequency responses. 

Identification of the second-order kernel due to 
pitching motions about the quarter-chord location was 
then performed using the same unit angle-of-ailack rate 
pulse input that was used for identification of the linear 
kernel. Nonlinear (transonic) aerodynamic unit pulse 
responses were computed at M=0.5,0.85, and 0.93 where 
strong shocks are present at the M=0.85 and the M=0.93 
conditions. Three terms of the second-order kernels for 
the three Mach numbers were then identified and 
discussed. The results are encouraging, although 
additional effort is required for determining the number of 
terms of the second-order kernel that are required for 
accurate prediction of nonlinear responses to arbitrary 
inputs and in verifying that the assumption of a weakly 
nonlinear second-order system is accurate for transonic 
aerodynamic responses. 
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Figure 1 The response of the linear system (N=0) 
and the nonlinear (bilinear) system due to the input 
shown in the inset for the example problem. 



Figure 2 Comparison of the analytical and numerically- 
identified linear unit impulse responses for the example 
problem. 



Figure 3 Numerically-identified second-order kernel 
at twenty values of time lag for the example problem. 



Figure 4 Comparison of analytical and numerically- 
identified second-order kernels at two values of time 
lag for the example problem. 



to pitch at the quarter-chord for the linear (flat plate) 
aerodynamic solution at M=0.85. 


C L and a, rad/sec 



Figure 6 Lift-coefficient response to the pitching motion 
about the quarter-chord location described by the angle-of- 
attack (inset) and angle-of-attack rate functions for the linear 
aerodynamic (fiat plate) solution at M=0.85. 
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Figure 7 Comparison of lift-coefficient responses, 
due to the exponential pitch pulse shown in fig. 6, 
for the CAP-TSD solution and the convolution of 
PI (fig. 5) and the angle-of-attack rate (fig.6). 



Figure 8 The nonlinear (thickness) lift-coefficient unit 
pulse responses PI, P2, PI 1, and P12 due to pitch at the 
quarter-chord location for the NACA64A010 rectangular 
wing at M=0.85 and y = 1 .4. 



Figure 9 Comparison of the lift-coefficient responses 
due to pitch at the quarter-chord location for the linear 
(fiat plate) case and the linear portion of the nonlinear 
(thickness) response case at M=0.85 and y = 1 .4. 



time steps 


Figure 10 Three terms of the second-order kernel of lift 
coefficient due to pitch about the quarter-chord location 
for the NACA64A010 rectangular wing at M=0.85 and 
y= 1.4. 



Figure 1 1 Comparison of nonlinear (thickness) lift- 
coefficient unit pulse responses due to pitch at the 
quarter-chord location for M=0.5, M=0.85, and M=0.93 
with y= 1.4. 



Figure 12 Three terms of the second-order kernel of lift 
coefficient due to pitch about the quarter-chord location 
for the NACA64A010 rectangular wing at M=0.50 and 
y= 1.4. 


12 






Figure 13 Three terms of the second-order kernel of lift 
coefficient due to pitch about the quarter-chord location 
for the NACA64A010 rectangular wing at M=0.93 and 
y= 1.4. 
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Figure 14 Relative nonlinearity for each of the 
three terms of the second-order kernel at M=0.5, 
M=0.85, and M=0.93. 



Figure 15 Comparison of linear (flat plate) and nonlinear 
(thickness) lift-coefficient responses due to the exponential 
pitch pulse shown in figure 6 at M=0.85 and y = 1 .4 (for 
nonlinear case). 
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